rm(list = ls())

library(tidyverse)
library(lfe)
library(broom)
library(stargazer)
library(janitor)

## Get data

df <- read_rds("data/Data_Clean.rds") %>% 
  mutate(year = as.numeric(year)) %>% 
  mutate(unit = coalesce(ags, vb_key)) %>% 
  mutate(treat_weakly = ifelse(treat_categ == 'schwach', 1, 0),
         treat_heavy = ifelse(treat_categ %in% c('schwer', 'sehr schwer'), 1, 0)) %>% 
  mutate(treat_weakly_post = treat_weakly * post,
         treat_heavy_post = treat_heavy * post) %>%
  filter(year %in% c('2017', '2021')) %>%
  filter(!is.na(unit)) %>% 
  mutate(treat_either = ifelse(treat_weakly == 1 | treat_heavy == 1,
                               1, 0)) 

## Define incumbent vote share 

df <- df %>% 
  mutate(incumbent_vs = case_when(land == 'NRW' ~ cdu_csu,
                                  land == 'RP' ~ spd),
         incumbent_vs = incumbent_vs*100) ## in p.p.  

## Difference the DV 

diff_df <- df %>%
  dplyr::select(one_of('year', 'unit', 'incumbent_vs', 'land', 'treat_weakly', 
                       'treat_heavy', 'treat_wounded_or_dead', 'county_id', 
                       'treat_either')) %>%
  pivot_wider(id_cols = c('unit', 'land', 'treat_weakly', 'treat_heavy', 
                          'treat_wounded_or_dead', 'county_id', 'treat_either'),
              values_from = 'incumbent_vs', 
              names_from = 'year') %>%
  clean_names() %>%
  mutate(diff_incumbent_vs = x2021 - x2017)

## Function to do this for the covars as well 

get_diff_var <- function(v) {
  
  out <- df %>%
    dplyr::select(one_of(c('unit', 'year', v))) %>%
    pivot_wider(id_cols = c('unit'),
                values_from = v, 
                names_from = 'year') %>%
    clean_names() %>%
    mutate(diff = x2021 - x2017)
  colnames(out)[4] <- v
  out %>% 
    dplyr::select(1, 4)
  
}

cov_list <- c("county_age", "county_unem", 
              "county_pop_tot", 
              "county_foreign_share", 
              "county_abitur_share")

cov_df <- cov_list %>% 
  lapply(get_diff_var) %>% reduce(left_join)

## Merge

diff_df <- diff_df %>% 
  left_join(cov_df)

## Table 3: estimation 

m1 <- felm(diff_incumbent_vs ~ treat_weakly + treat_heavy | land | 0| county_id,
           data = diff_df)

m2 <- felm(diff_incumbent_vs ~ treat_wounded_or_dead | land | 0| county_id,
           data = diff_df)

m3 <- felm(diff_incumbent_vs ~ treat_weakly + treat_heavy | 0| 0| county_id,
           data = diff_df %>% 
             filter(land == 'RP'))

m4 <- felm(diff_incumbent_vs ~ treat_wounded_or_dead | 0| 0| county_id,
           data = diff_df %>% 
             filter(land == 'RP'))

m5 <- felm(diff_incumbent_vs ~ treat_weakly + treat_heavy | 0 | 0| county_id,
           data = diff_df%>% 
             filter(land == 'NRW'))

m6 <- felm(diff_incumbent_vs ~ treat_wounded_or_dead | 0| 0| county_id,
           data = diff_df%>% 
             filter(land == 'NRW')) 

mlist <- list(m1,m2,m3,m4,m5,m6)

stargazer::stargazer(list(m1,m2,m3,m4,m5,m6), keep = 'treat',
                     covariate.labels = c('Flooding: weakly affected',
                                          'Flooding: highly affected',
                                          'Any wounded or dead'),
                     style = 'ajps',
                     dep.var.labels = 'DV: Change in incumbent vote share (p.p.)',
                     model.numbers = F, 
                     keep.stat = c('n', 'rsq'))

## REPEAT FOR TABLE A.7 (including covariates)

## Estimate 

m1 <- felm(diff_incumbent_vs ~ treat_weakly + treat_heavy +
             county_age + county_unem +
             county_pop_tot +
             county_foreign_share +
             county_abitur_share| land | 0| county_id,
           data = diff_df)

m2 <- felm(diff_incumbent_vs ~ treat_wounded_or_dead+
             county_age + county_unem +
             county_pop_tot +
             county_foreign_share +
             county_abitur_share | land | 0| county_id,
           data = diff_df)

m3 <- felm(diff_incumbent_vs ~ treat_weakly + treat_heavy+
             county_age + county_unem +
             county_pop_tot +
             county_foreign_share +
             county_abitur_share | 0| 0| county_id,
           data = diff_df %>% 
             filter(land == 'RP'))

m4 <- felm(diff_incumbent_vs ~ treat_wounded_or_dead +
             county_age + county_unem +
             county_pop_tot +
             county_foreign_share +
             county_abitur_share| 0| 0| county_id,
           data = diff_df %>% 
             filter(land == 'RP'))

m5 <- felm(diff_incumbent_vs ~ treat_weakly + treat_heavy +
             county_age + county_unem +
             county_pop_tot +
             county_foreign_share +
             county_abitur_share| 0 | 0| county_id,
           data = diff_df%>% 
             filter(land == 'NRW'))

m6 <- felm(diff_incumbent_vs ~ treat_wounded_or_dead +
             county_age + county_unem +
             county_pop_tot +
             county_foreign_share +
             county_abitur_share| 0| 0| county_id,
           data = diff_df%>% 
             filter(land == 'NRW')) 

mlist <- list(m1,m2,m3,m4,m5,m6)

stargazer::stargazer(list(m1,m2,m3,m4,m5,m6), keep = 'treat',
                     covariate.labels = c('Flooding: weakly affected',
                                          'Flooding: highly affected',
                                          'Any wounded or dead'),
                     style = 'ajps',
                     dep.var.labels = 'DV: Change in incumbent vote share (p.p.)',
                     model.numbers = F, 
                     keep.stat = c('n', 'rsq'))




